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ABSTRACT 

For a simple Jeffcott model with deadband, we develop a method; to 
determine stability margins by an analysis of the Discrete Fourier transform 
of the system response. The viability of this method is demonstrated by 
means of numerical simulations. We also expktin the circular behavior of the 
system response on the stability boundary^ 
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INTRODUCTION 


In this report we develop a method to determine stability margins for a 
simple Jeffcott model with deadband. We consider an unbalanced, uniform, 
flexible shaft supported by bearing, that rotates along the x axis. If m is 
the mass of the shaft, y and z describe the displacement of the center of the 
shaft, r = (y 2 + z 2 y^ 2 ,h(t) = 1 if r < S and h(t ) = S/r if r > 6, then our 
model is represented by the following system of coupled nonlinear differential 
equations: 


y" + Cy' + [A + K( 1 - h{t))]y + Bz = /,(() . , 

z" + Cz' -By + [A + K( 1 - A(i)))z = h(t), ' ' 

where the coefficients have the following physical interpretation: C = C a /m, A 
K a /m,K = Kb/m,B = Q a /m where C a is the seal damping, K s the seal 
stiffness, Kb the bearing stiffness, and Q a denotes the cross-coupling stiff- 
ness of the seals. We assume that fi(t) and f 2 {t) are arbitrary continu- 
ous and bounded functions. In the particular case in which the shaft ro- 
tates with angular velocity w, and no sideloads or other factors are present, 
fi(t ) = d m w 2 cos wt,f 2 (t) = d m w 2 siniof, and d m is the displacement of the 
shaft’s center of mass from the geometric center. A sketch of our model can 
be seen in Fig. 1, with the caveat that the seals are not shown. To say that 
our model is stable, is the same as saying that all the solutions of the system 
of differential equations (1) are bounded. 

THEORETICAL ANALYSIS 

This is a review of results obtained in [1], to which the reader is referred for 
additional details. This paper was motivated by previous work of Day [2], 
Yamamoto [3], and Childs [4], See also Rowan [5]. 

Setting M = A + K — iB, n = y + iz,q = hfi, f = f\ + if 2 , multiplying 
the second equation in (1) by i and adding it to the first one, we see that (1) 
is equivalent to 


fi" + <V + Mu - Kq = f{t). (2) 

In order to study the stability of (2) we first study the stability of the lin- 
earized system 


(3) 


li' + Cp' + Mp = f(t). 

It can be shown that the solutions fit of (3) have the form 

fit = ci exp^t) + c 2 exp(A 2 t) + fi P , (4) 

where fi p is a particular solution of (3), Ai and A 2 are the roots of A 2 + 
C A + M = 0, and Ci,c 2 are arbitrary constants. It is also known that 
Aj = a+ift, A 2 = a'—ift, where a , a' and ft are real and a' < a. Clearly ft/2ir 
is the natural frequency of the linear system (3). The following inequalities 
will be useful in the sequel: 

If a < 0, then ^ < ft < \J A -f K , and a 1 < 0 
If a = 0, then ^ = ft = \J A -f K , and a' < 0 (5) 

If a > 0, then ® > ft > y/A + K 

We have shown in [1] that, if G(t) = (Ai — A 2 ) _1 [exp(Ai<) — exp(A 2 t)], then 
every solution fi of (2) can be represented in the form fi(t) = fii(t) + P(t), 
where 


P(t) = K f G(t — x ) • q(x)dx. 

Jo 

Since |g(<)| < 6, it readily follows that P(t) is bounded whenever a < 0. An 
inspection of (4) readily shows that this is the same condition that guaran- 
tees the boundedness of fit. In other words, (2) is stable if and only if the 
linearized system (3) is stable. However, since the coefficients of (2) are in 
general not very well known, this conclusion does not have a straightforward 
application. We shall now see that, nevertheless, it can be used in an analysis 
of the discrete Fourier transform of the system response fi(t). 

Let Ofc >n = nd + k = 0, . . . , m — 1, d = m£, m and £ constant, and 

m-1 

D n fi{ A) = m" 1 M«fc,n) ex P(-^ a it,n)- 

k=0 

Assume fi(t) = transients -f YikLo hc ,<rkt , with YlkLo IM < c* 0 , and let A j,j = 
0, . . . , N be given. Then, for any e > 0 there are £ and m such that 

= -tW-V) + KU + K- A?) + 
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and | Ej,n \ < £,j = 0, ...,N,n = 0,1,2,... We have therefore obtained the 
approximate identity 


DM A,) = D^XJ + KW+K-W+iiCXj-By-'DM^i = 0 N 

(6) 

Note that we have not assumed that the signal is band limited. Most of the 
remainder of this report deals with material that will appear in [6]. 

APPLICATIONS 

We define the (theoretical) stability margin SMT of the system represented 
by (2), by means of the identity 

SMT = |(A + K - ft 2 ) + i(Cft - B)\/K 

In view of (5), SMT — ► 0 as a — ► 0. Thus, if D n g, p (ft) = 0 (which would 
happen if for instance f(t ) = A 0 exp(iwt) and ft ^ w ), we infer from (6) that 
SMT could be computed as the ratio of \D n q(ft)\ and \D n g(ft)\. 

In practice, the value of ft is of course unknown. However, if D n q(ft) does 
not vanish as a — ► 0, the right hand member of (5), evaluated at A_, = ft, 
would diverge as a — ► 0. This would mean that the graph of |Z) n /i(A)j would 
exhibit a spike at some frequency f n t (which we shall call the nonlinear 
natural frequency) that approaches /?/27r as a — ► 0 (i.e. the linear and 
nonlinear natural frequencies coalesce) and that the magnitude of the spike 
would generally increase as a decreases. Thus, the ratio \D n q(\) / D n y,(X)\ 
evaluated at A = 2irf n ( is an approximation to SMT and should yield an 
estimate of the stability margin of the system. 

EXAMPLES 

We consider (2) with initial conditions /z(0) = £,//(0) = 0, and A = 0, C = 
240,6 = 0.0000285, K = 1,305, 000, d m = 0.000060915, w = 1,000 irs,B = 
60 to, where s varies. The linear natural frequency ft/2 t will be denoted by / 0 . 
The forcing frequency /i clearly equals w/2rc, or 500s. The critical frequency 
f c is the value of / 0 whe n o = 0. From (4) we have f c = \f A + K /( 2n) = 
181.8 Hz. Since ® = y/A + K, solving for s we readily deduce that this 
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frequency corresponds to s c = 1.4545. This is the same example that we 
studied in [1]. 

In Table 1 we compare linear and nonlinear natural frequencies for various 
values of s, and different sampling intervals and sample sizes, and verify the 
assertion that both frequencies coalesce. We also observe that the location of 
the nonlinear frequency is practically independent of sampling interval and 
sample size. 

In Table 2 we estimate stability margins, defined as the ratio \D n q/D n pi\ 
computed at the location of the nonlinear natural frequency. We note that 
for each sampling interval and sample size the computed stability margins 
decrease as a increases. Tables 1 and 2 have been obtained from numerical 
solutions of the examples using a fourth order Runge-Kutta method followed 
by a Fast Fourier Transform algorithm. This procedure simulates an analy- 
sis of samples of the system response obtained through displacement probes. 
However, in many situations these probes are unavailable or yield inaccurate 
measurements and therefore accelerometer data becomes the only reliable 
kind available. To obtain the Discrete Fourier transform of (i from the Dis- 
crete Fourier transform of pi", we can use the approximate identity 


D„ii"( A) = -A’A^A), (7) 

which can be readily deduced from the mathematical analysis developed in 

[1]. 

In Fig. 2(a) we plot the absolute value of D n u for s = 1.3 and 1001 points 
in the interval 2.56 < t < 3.56. The same plot, obtained from D n pi"( X) using 
(7), is shown in Fig. 2(b). Note that both plots are identical, except at low 
frequencies. Fig. 3 shows a plot of the absolute value of D n q (obtained by a 
simulation of a displacement probe), for the same range of parameters that 
was used in Fig. 2. 

Formula (7) suggests a method to determine stability margins from ac- 
celerometers: Given a sample of the acceleration pi" of the system response, 
compute D n pt(kj),j = 0, . . . , N using (7). To find D n q , do an inverse Fourier 
transform to estimate From pi it is then possible to estimate q(t). For 
each value of t , this estimate will be very inaccurate. However, it is rea- 
sonable to expect that for sufficiently large N these errors will even out, 
yielding a reasonable good estimate of D n q. The feasibility of this method 
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is demonstrated in Table 3. 


ANNULAR VS. CIRCULAR BEHAVIOR 

It was observed by Day [2] that near the stability boundary the system 
response has circular behavior, and there has been some intertest in finding 
an explanation for this phenomenon. The explanation is, in fact, very simple: 
Fig. 4 shows a plot of \D n fi\ on the stability boundary. Clearly the nonlinear 
natural frequency f 0 dominates all other frequencies, and therefore D n /i( A) = 
A 0 «5(A — / 0 ), where 8 denotes Dirac’s delta function. Inverting, we have 
Dnfi(X) = Aiexp(if 0 t). Fig. 5 exhibits the circular behavior of the system 
response on the stability boundary. 

MAJOR POINTS 

* Proof that nonlinear frequency approaches the frequency of instability 
at the stability boundary 

* PSD of displacement can be obtained from PSD of acceleration 

** Given knowledge of the location of the nonlinear spike and of the psd 
of acceleration, a real time direct calculation of stability margin can 
be made 

* An explanation of circular solutions vs. annular 

* Methods have been proven only on a simple, low order system 

* Physical proof of the method on a real system is the next logical step 
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TABLES AND FIGURES 



s 

alpha 

fO 

d=l 

d=128 

1.3 

-12.63 

181.61 

163 

162.74 

1.35 

-8.54 

181.68 

169 

168.78 

1.4 

-4.45 

181.74 

175 

175.00 

1.45 

-0.37 

181.81 

181 

181.25 

1.4545 

-4.51E-4 

181.81 

181 

181.80 

1.5 

3.71 

181.88 

182 

181.51 


Table 1. Comparison of Linear and nonlinear natural frequencies. Column 
three shows linear natural frequencies. Columns four and five show observed 
nonlinear natural frequencies for intervals of the form (2.56,2.56 + d). The 
results in column four were obtained using both 2 U and 2 18 points. The 
values obtained were the same. For column five, 2 18 points were used. 
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s 

alpha 

1.3 

-12.63 

1.35 

-8.54 

1.4 

-4.45 

1.45 

-0.37 

1.4545 

-4.51E-4 

1.5 

3.71 

1.6 

11.85 


n = 2 11 


SMT 

d=l 

2.22E-2 

1.99E-1 

1.50E-2 

1.38E-1 

7.83E-3 

7.35E-2 

6.47E-4 

8.71E-3 

7.93E-7 

5.01E-3 

6.54E-3 

2.06E-7 

2.09E-2 

2.48E-19 


n = 2 18 

n = 2 18 

d=l 

D- 

II 

h-* 

bO 

00 

1.99E-1 

1.99E-1 

1.38E-1 

1.38E-1 

7.35E-2 

7.36E-2 

8.71E-3 

6.18E-3 

4.99E-3 

1.75E-4 

2.06E-7 

2.07E-212 

2.49E-19 



Table 2. Determination of stability margins. The third column lists 
the theoretical stability margins. Columns four, five and six list stability 
margins computed from intervals of the form [2.56,2.56 + d], with sample 
size n, sample rate 1/n, and d and n as indicated. 






d = 1 

d= 1 

d = 128 

s 

a 

SMT 

n = 2 12 

n = 2 18 

n = 2 18 

1.3 

-12.63 

2.22E-2 

4.84E-3 

1.92E-3 

2.03E-1 

1.35 

-8.54 

1.50E-2 

2.05E-3 

8.92E-4 

5.21E-2 

1.4 

-4.45 

7.83E-3 

3.80E-2 

1.59E-2 

6.78E-2 

1.45 

-0.37 

6.47E-4 

4.57E-4 

4.15E-4 

5.08E-3 

1.4545 

-4.51E-4 

7.93E-7 

4.25E-5 

1.17E-5 

2.82E-5 


Table 3. Determination of stability margins from simulated accelerom- 
eter data. Columns four, five and six list stability margins computed from 
intervals of the form [2.56, 2.56 +d], with sample size n, sample rate 1/n, 
and d and n as indicated. 
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FIG. 1 
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5 = 1 . 3 , 2.56 <*< 3 . 56 , 1001 ^ 5 . 

Fig. 2(a) 
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s = 1.3, 2.56 < t < 3.56 

PSD plot of y obtained from PSD plot of y" 

Fig. 2(b) 
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PSD plot for q, s = 1.3 

2,56 < t < 3.56, 1001 points 

Fig. 3 
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PSD plot for s = 1.454505522, 2.56 < t < 3.56 
value at f Q : 4.55*10~"* 

value at f^: 5.9*10 - ^ 

Fig. 4 
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TOO TOO- 


o 


- 0.01 0.01 

2.56 < t < 2.81, 2^ points, 

s = 1.454505522, a = -2,59*10 -8 

min | p(t) | = 4.78*10 -3 , max|y(t)| = 5.33*10 _3 

Fig. 5 
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